
res<-clinical_univ_res
prefix<-"uni_clinical_info"
forest_plot <- function(res,prefix){
  result <- cbind(rownames(res),res)
  colnames(result)[1] <- "features"
  write.table(result,paste0(opt$outdir,'/',prefix,'_cox_univ_res.txt'),sep="\t",quote=F,row.names=F)
  res$univ_beta <- round(res$univ_beta,3)
  res$univ_HR <- round(res$univ_HR,3)
  res$`univ_95% CI for HR` <- sapply(res$`univ_95% CI for HR`,function(x){paste0(signif(as.numeric(unlist(strsplit(x,'-'))),2),collapse = "-")})
  res$univ_p.value <- round(res$univ_p.value,3)
  res$mean <- as.numeric(as.character(res$univ_HR))
  res$lower <- as.numeric(as.character(sapply(res$`univ_95% CI for HR`,function(x){unlist(strsplit(as.character(x),'-'))[1]})))
  res$upper <- as.numeric(as.character(sapply(res$`univ_95% CI for HR`,function(x){unlist(strsplit(as.character(x),'-'))[2]})))
  text <- rbind(colnames(res)[1:4],as.matrix(res[1:4]))
  rownames(text)[1] <- 'features'
  text <- cbind(rownames(text),text)
  plotHR(text,res,prefix)
}

plotHR <- function(text,res,prefix){
  text[1,] <- gsub("univ_|multi_","",colnames(text))
  pdf(paste0(opt$outdir,'/',prefix,'_forestplot.pdf'),width=10, height=5)
  # pdf("forestplot.pdf",width=960, height=(1+dim(res)[1])*50)
  colinfo <- lapply(1:nrow(res),function(x){ifelse(x %% 2==1,
                                                   return(gpar(col="gray",lineend="butt",columns=c(2:6))),return(gpar(col='gray',lineend="butt",columns=c(2:6))))})
  
  names(colinfo) <- 1:nrow(res)+1
  forestplot(as.matrix(text),mean=c(NA,log2(res$mean)),lower = c(NA,log2(res$lower)),upper = c(NA,log2(res$upper)),
             graph.pos = 5,
             graphwidth = unit(60,'mm'),
             fn.ci_norm="fpDrawDiamondCI",
             lineheight = unit(12,'mm'),
             line.margin = unit(5,'mm'),
             colgap = unit(4,'mm'),
             zero = 0,cex=0.9,
             col=fpColors(box="#DC143C", lines="#FFC0CB", zero = "black"),
             txt_gp = fpTxtGp(label=gpar(cex=1.25),
             ticks=gpar(cex=1.1),xlab=gpar(cex=1.2),title = gpar(cex=1.2)),
             lwd.ci=1.5,ci.vertices = T,title="log2 Hazard Ratio",boxsize = 0.4,hrzl_lines = colinfo)
  dev.off()
}

res<-clinical_multiv_res
prefix<-"multi_clinical_info"

forest_plot_multi <- function(res,prefix){
  result <- cbind(rownames(res),res)
  colnames(result)[1] <- "features"
  write.table(result,paste0(opt$outdir,'/',prefix,'_cox_multi_res.txt'),sep="\t",quote=F,row.names=F)
  res$multi_beta <- round(res$multi_beta,3)
  res$multi_HR <- round(res$multi_HR,3)
  res$`multi_95% CI for HR` <- sapply(res$`multi_95% CI for HR`,function(x){paste0(signif(as.numeric(unlist(strsplit(x,'-'))),2),collapse = "-")})
  res$multi_p.value <- round(res$multi_p.value,3)  
  res$mean <- as.numeric(as.character(res$multi_HR))
  res$lower <- as.numeric(as.character(sapply(res$`multi_95% CI for HR`,function(x){unlist(strsplit(as.character(x),'-'))[1]})))
  res$upper <- as.numeric(as.character(sapply(res$`multi_95% CI for HR`,function(x){unlist(strsplit(as.character(x),'-'))[2]})))
  text <- rbind(colnames(res)[1:4],as.matrix(res[1:4]))
  rownames(text)[1] <- 'features'
  text <- cbind(rownames(text),text)
  plotHR(text,res,prefix)
}


plotHR <- function(text,res,prefix){
  text[1,] <- gsub("univ_|multi_","",colnames(text))
  pdf(paste0(opt$outdir,'/',prefix,'_forestplot.pdf'),width=10, height=5)
  # pdf("forestplot.pdf",width=960, height=(1+dim(res)[1])*50)
  colinfo <- lapply(1:nrow(res),function(x){ifelse(x %% 2==1,
                                                   return(gpar(col="gray",lineend="butt",columns=c(2:6))),return(gpar(col='gray',lineend="butt",columns=c(2:6))))})
  names(colinfo) <- 1:nrow(res)+1
  forestplot(as.matrix(text),mean=c(NA,log2(res$mean)),lower = c(NA,log2(res$lower)),upper = c(NA,log2(res$upper)),
             graph.pos = 5,
             graphwidth = unit(60,'mm'),
             fn.ci_norm="fpDrawDiamondCI",
             lineheight = unit(12,'mm'),
             line.margin = unit(5,'mm'),
             colgap = unit(4,'mm'),
             zero = 0,cex=0.9,
             col=fpColors(box="#DC143C", lines="#FFC0CB", zero = "black"),
             txt_gp = fpTxtGp(label=gpar(cex=1.25),
                              ticks=gpar(cex=1.1),xlab=gpar(cex=1.2),title = gpar(cex=1.2)),
             lwd.ci=1.5,ci.vertices = T,title="log2 Hazard Ratio",boxsize = 0.3,hrzl_lines = colinfo)
  dev.off()
}
